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Abstract 



(N 
> 

We investigate {d + l)-dimensional quasilinear systems which are integrable by 
the method of hydrodynamic reductions. In the case d> 3 we formulate a conjec- 
^ ture that any such system with an irreducible dispersion relation must be linearly 

degenerate. We prove this conjecture in the 2-component case, providing a complete 
classification of multi-dimensional integrable systems in question. In particular, our 
O results imply the non-existence of 2-component integrable systems of hydrodynamic 

^ type for d>6. 

IT^ In the second half of the paper we discuss a numerical and analytical evidence 

j_] for the impossibility of the breakdown of smooth initial data for linearly degenerate 



systems in 2 -|- 1 dimensions. 
MSG: 35L40, 35L65, 37K10. 
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1 Introduction 



In this paper we discuss {d + l)-dimensional quasilinear systems of the form 

d 

+ ^Ai(u)u^, = 0, (1) 

i=l 

where t,a;* are independent variables, u is an ra-component column vector and Aj(u) are 
n X n matrices. Systems of this type play an important role in continuum mechanics 
and physics, and there is a vast literature covering their analytical, geometric and applied 
aspects, see e.g. [71 EU [36l ISTJ [191 El US]- In this paper we assume that the dispersion 
relation det(/i/ + a*Aj) = defines an irreducible algebraic variety (for n = 2, an irre- 
ducible quadric) in the variables /i, a*. Our prime goal is the classification of integrable 
systems of the form ([T]) in any dimension. 

For d = 1, any two-component system of the form 

ut + A{u)u^ = 

is linearizable by a hodograph transformation which interchanges the dependent and in- 
dependent variables and, therefore, is automatically integrable. In particular, any such 
system possesses an infinity of conservation laws and commuting flows. 
For d = 2, 

Ut + A{u)u^ + B{u)uy = 0, (2) 

the situation is more subtle. Even the very concept of integrability for multi-dimensional 
systems of the form ([T]) requires an alternative approach. Such approach, based on the 
method of hydrodynamic reductions, was outlined in [13]. It was proposed to call a 
(2 -|- l)-dimensional system integrable if it possesses infinitely many multi-phase solutions 
of the form u = u{R^, R^) where the phases R^{x,y,t) satisfy a pair of commuting 
hydrodynamic type systems 

R\ = X'{R) Rl, Ry = Rl; (3) 

we emphasize that the number of phases N is allowed to be arbitrary. In what follows, 
we will refer to the equations ^ as A^-component hydrodynamic reductions. Multi-phase 
solutions of this type, known in gas dynamics as nonlinear interactions of planar simple 
waves, date back to [HI [6l |3ll [20]. Their first appearance in the context of integrable 
systems and, in particular, the dispersionless Kadomtsev-Petviashvili hierarchy, is due to 
[TTJ [TS]. The requirement of the existence of multi-phase solutions imposes strong con- 
straints on the matrices A{u) and -B(u). In the two-component case these integrability 
conditions were derived in [H], and are included in Sect. 4.1 of this paper. In partic- 
ular, they imply that any integrable system of this kind is necessarily conservative, and 
possesses a dispersionless Lax pair. It was demonstrated in [301 131] that the integrability 
conditions can be resolved in generalized hypergeometric functions. 

Remarkably, for d > 3 the situation simplifies dramatically: the requirement of the 
integrability proves to be too strong, so that very few examples survive in many dimen- 
sions. According to the definition proposed in [15], a {d + l)-dimensional system of the 
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form ([T]) is said to be integrable if it possesses infinitely many A^-component hydrody- 
namic reductions parametrized by {d — 1)N arbitrary functions of one variable. Here 
by a hydrodynamic reduction we mean a decomposition of a (c? + l)-dimensional system 
into d commuting systems of the form (|3|. Although the derivation of the integrability 
conditions based on the method of hydrodynamic reductions leads to a quite complicated 
analysis, there exists a simple way to bypass lengthy calculations. It is based on the 
requirement that all traveling wave reductions of a (rf + l)-dimensional integrable system 
to (2 + l)-dimensional systems should be themselves integrable. Since the integrability 
conditions in 2 + 1 dimensions are explicitly known, this provides a set of strong con- 
straints which are therefore necessary for the integrability. In fact, in our case they turn 
out to be sufficient. Our main result is the following 

Theorem 1 There exist no two-component integrable systems of the form ^ in dimen- 
sions 6 + 1 and higher, /n 5 + 1 dimensions, there exists a unique integrable system. Up 
to linear transformations of the independent variables, it can be brought to the form 

Vt = + WVr — VWr, Wz = Vy + VWg — WVg. (4) 

/n 4 + 1 and 3 + 1 dimensions, any two-component integrable system can be obtained as 
a traveling wave reduction of the system 

The proof of this result is given in Sect. 4.2. The system Q first appeared in p5], see 
Sect. 2 for more details. An important property of this example is its linear degeneracy. 
Recall that a matrix A{u) is said to be linearly degenerate if its eigenvalues (assumed real 
and distinct) are constant in the direction of the corresponding eigenvectors. Explicitly, 
Lj-iX^ = 0, no summation, where L^i is the Lie derivative in the direction of the eigenvector 
r*, and Ar^ = AV*. A system of the form ([l]) is said to be linearly degenerate (totally 
linearly degenerate in the terminology of [2lj), if any matrix A{u) = a^Ai{u) is linearly 
degenerate for any values of the constants a*. Theorem 1 and other existing examples 
support the following conjecture. 

Conjecture. For d > 3, any n-component (c? + 1)- dimensional integrable system of the 
form ([7]j with an irreducible dispersion relation must be linearly degenerate. 

According to the conjecture of Majda [21], p. 89, linearly degenerate systems are quite 
exceptional from the point of view of the global solvability of the Cauchy problem: the 
shock wave formation, which is a standard scenario for genuinely nonlinear systems, never 
happens for smooth initial data. Although this fact is well-established in 1 + 1 dimen- 
sions ESI ESI EZ] , the mult i- dimensional situation is still out of reach. We refer to the 
recent work [251 EHl ET] which utilizes a novel version of the inverse scattering transform 
to investigate the breakdown phenomena for some (2 + l)-dimensional dispersionless in- 
tegrable models. In Sect. 3 we provide some further analytical and numerical evidence in 
support of Majda's conjecture. As an example, we choose a particular (2 + l)-dimensional 
traveling wave reduction of the system (|4]), namely, 

vt + w^ = 0, Wt + Wy + vw^ - wvx = 0, (5) 

see Sect. 3.1 for the discussion of this example. In a somewhat different form, it appeared 
in [39j. In Sections 3.2-3.3, applying the method of hydrodynamic reductions and other 
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simple techniques, we construct some global exact solutions of the system (|5]). In Sect. 3.4 
we summarize the results of numerical simulations based on the Friedrichs-symmetrized 
representation of the system ([s]). These results clearly indicate that hump-like initial data 
tend to spread out, and do not break down in finite time. 



2 Examples of mult i- dimensional integrable systems 

In this section we list some examples of integrable l)-dimensional quasilinear systems 
of the form ([T]) for d > 3. All such systems turn out to be linearly degenerate, namely, an 
arbitrary matrix of the linear pencil ^(u) = aMj(u), where a* are arbitrary constants, 
must be linearly degenerate (totally linearly degenerate systems in the terminology of [24]). 
We recall that there exists a useful criterion of linear degeneracy which can be summarized 
as follows. Given a matrix A{u), let us introduce its characteristic polynomial, 

det{XE - A{u)) = A" + /i(u)A"-i + /2(u)A"-2 + . . . + /„(u). 

The condition of linear degeneracy is given by [12] , 

V/i + V/2 + . . . + v/„ = 0, (6) 

where V is the operator of the gradient, V/ = {§^, • • • , -§^), and A'' denotes fc-th power 
of the matrix A. Our first example plays a key role in the classification of two-component 
multi-dimensional quasilinear systems. 

Example 1. The six- dimensional two-component integrable system Q, 

Vt = Wx + WVr — VWr, Wz = Vy + VWs — WVs, 



first appeared in . It can also be obtained as the condition of commutativity of two 
parameter-dependent vector fields, 

[dz — vds — Xdx + Xvdr, dy — wds — Xdt + Xwdr] = 0. 

It was demonstrated in |15j that system Q possesses infinitely many n-component hydro- 
dynamic reductions parametrized by 4n arbitrary functions of one variable. This example 
is a result of the following construction. In [14J, we proposed a characterization of two- 
component (2 -|- l)-dimensional integrable systems of the form 

The integrability conditions constitute a complicated overdetermined system of second 
order PDEs for the 2x2 matrices A and B which are presented in Sect. 4.1. In the 
particular case when A is taken to be linearly degenerate. 



A=\ \, (7) 




the integrability conditions imply the following generic expression for B, 

/ ^-^-^ — aw 

111 — 01 



B 



f(w) 2 fiv) 

w—v w~v 



^ v—w v~w 

3 I /3.,2 



where / is a cubic polynomial, f{v) = av + (3v + 'jv + S, and a,(3,'j,5 are arbitrary 
constants. A remarkable property of this example is that any matrix in the linear pencil 
B + fiA is also linearly degenerate. Explicitly, one has 

B = 6Bi+ + /3B3 + aBi = 

( 1 1 \ / W V \ 

w—v w—v I I /-,/ I w—v w—v \ _i_ IQ I w—v w—v 1 I „, 
V-W V-W J \ V-W V-W J \ v-u) v-w I \ v-w v-w 

Let us introduce the (5 + 1) -dimensional system 

Notice that an arbitrary linear combination of matrices A and i?2, -B3, -B4 is automati- 
cally linearly degenerate. This can be verified directly using the criterion (|6|. In the new 
variables V = v \ w, W = vw, our system reduces to 

Vt + W^ + WVr - VWr = 0, Wt + VW^ - WV^ + Vy + W, + VWs - WVs = 0, 

taking a fully symmetric form Q after an obvious linear transformation of the independent 
variables. We would like to emphasize that there was no a priori reason for the (5 + 1)- 
dimensional system obtained in this way to be integrable: this is a miracle which is due 
to the linear degeneracy. Notice that a particular dimensional reduction of the system 

mt = n^, nt = my + mn^ - nm^, (9) 

known as the Pavlov equation, has been extensively discussed in the context of hydrody- 
namic chains and the 'universal hierarchy' [331 ESI SOl l27] - 

Example 2. The following four-dimensional three-component integrable system appeared 
in ^llj in the context of hyper-Hermitian conformal structures, 

AN, + 5iNyi + B2^y2 + fi3Nj,3 = 0, 

here N is a three-component column vector, N = {Ni, N2, N^Y, and A, Bi, B2, B3 are 
3x3 matrices, 

/ 1 -N3 N2 \ /ooo\ „ /oo-i\ „ /oio 

A = { N:i I -Ni ] . = 1 , i?2 = , i?3 = -1 

\-N2 1/' Vo-lOy' VlOO/' VOOO 

Multiplying by A~^ one can bring this system to an evolutionary form, moreover, any 
linear combination of matrices A~^Bi is linearly degenerate. Notice, however, that the 
dispersion relation of this system is not irreducible. 

A general scheme for constructing multi-dimensional integrable systems based on the 
commutativity of vector fields was discussed in [3]. A broad class of new examples of 
integrable linearly degenerate systems in 2 + 1 dimensions was constructed in 
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3 Global solutions of linearly degenerate systems 



In 1 + 1 dimensions, linearly degenerate systems are known to be quite exceptional from 
the point of view of the solvability of the Cauchy problem: generic smooth initial data 
do not develop shocks in finite time [351 EHl [231 EZ]- The conjecture of Majda [21], p. 
89, suggests that the same statement should be true in higher dimensions, namely, for 
linearly degenerate systems the shock formation never happens for smooth initial data. 
For a particular case of the Pavlov equation, this conjecture was confirmed in [27] based 
on the novel version of the inverse scattering transform. In this section we discuss another 
dimensional reduction of the system Q, namely, the system (|5|, 

Vt + Wa: = 0, Wt + Wy + VWx - WVx = 0. 



In contrast to the Pavlov equation (|9j), this system has a well-posed initial value problem. 
We also refer to [H |5] for a detailed discussion of the 10 x 10 linearly degenerate Born- 
Infeld system which is a nonlinear version of Maxwell's equations. We hope that our 
linearly degenerate example will be more appropriate for establishing the global existence 
type results. 

Some useful properties of system ^ are discussed in Sect. 3.1: these include the 
conservative representation and the Friedrichs-symmetrized form, which is most suitable 
for numerical simulations. In Sect. 3.2 and 3.3 we construct globally defined exact 
solutions of system ([s]) by reducing it to known integrable equations in 1 + 1 dimensions. 
This includes reductions to the sinh-Gordon and KdV equations (Sect. 3.2), as well as 
the reduction to a pair of A^-component linearly degenerate systems (Sect. 3.3). Sect. 3.4 
contains numerical simulations of the system ([s]) which support the conjecture that smooth 
initial data do not break down in finite time. We emphasize that all these properties are 
not specific to the system ([s]), and hold for any linearly degenerate integrable system. 



3.1 Properties of System ([5j) 

One can show that system ([s]) possesses exactly three zero order conservation laws, 

Vt + W:c = 0, 
(v"^ - W)t + {VW)^: -Wy = Q, 

It can be written in a Friedrichs-symmetrized form, 
-w \ / f \ / —w \ I V 



1 / \ w / ^ \ —w V 




0, (10) 



which shows that one should expect at least finite time existence results in the region 
w < 0. Introducing the variables p = —2v/w'^, q = one can rewrite in the 

Godunov form [T^ . 

{fp)t + {9p)x + {hp)y = 0, {fq)t + {gq)x + {hq)y = 0, 
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where 

Remark. The system ^ possesses a Lax pair which can be represented as the condition 
of commutativity of two vector fields, 

[dt + dy + {v + X)d,, wd, - Xdt] = 0. (11) 

Lax pairs of this kind form a basis of the novel generalization of the inverse scattering 
transform [251 ESI [2Z] • 



3.2 Reductions to integrable equations 

The system ([s]) possesses a number of remarkable reductions to soliton equations in 1 + 1 
dimensions. Below we describe two particular reductions, to the sinh-Gordon and the 
KdV equations, respectively. For convenience, we make a change of variables x,y,t 
x,y,T = t — y which brings system ([S]) into the form 

Vr + Wx = 0, Wy + VWx — WVx = 0. (12) 



Reduction to the sinh-Gordon equation. Applying the procedure of to the Lax 
pair (11) one can obtain a class of exact solutions of the system (12) in the form 

V = Uy cos X, ti; = — sinh u sin a; + cosh M, (13) 

where u{y, r) satisfies the sinh-Gordon equation 

Uyr = sinh u. 

In the original variables, it takes the form Uyt + Uu = sinhu. Taking a traveling wave 
solution in the form 



u 



arccosh (e — (e — 1) sn^[cr(y — ct) \ m]) sign ^cos '^^) ^ 



r = 2.-A-(m). --^^^y ™ = ^. (14) 

where sn[-|m] is a Jacobi elliptic function and K{m) is the complete elliptic integral of the 
first kind (e.g. [I]), < c < 1 and e > 1, one obtains bounded doubly periodic solutions 
v{x, y, t) = (uy + Ut) cos X and w(x, y,t) = — sinh -u sin x + cosh u of the system ^ shown 
in Fig. [ijfor c = 1/3 and e = 1.1 (plotted using Mathematica pH])- These solutions remain 
non-singular for all values of t. 

Reduction to the KdV equation. Let v{x,y,T) be a solution to a pair of commuting 
flows ^ 

Vt = (V5^)x, Vy = -Vxxx - 3VVx, 

7 



Figure 1: Solution (13) 



(14) to the system ([sj for v and at t = 0. 



where satisfies the system 



Notice that the last two equations constitute a Lax pair for the KdV equation (at zero 
value of the spectral parameter), and the consistency condition Vry = Vyr holds identically 
modulo the equations for ip. Setting 



(15) 



one can see that v{x,y,T) and w{x,y,T) solve the system (12). Exact non-singular multi- 
soliton solutions of these equations can be constructed in the form 



lnW{i!i, ...,V'n), f 



1 WjiP,, 1) 

l^(^i,...,V'n) ' 



(16) 



here W denotes the Wronskian of n functions ipi where 4'2j-i = cosh(fc2j_ia; + kl^j^iV + 
T /k2j-i) and ip2j = smh{k2jX + k^jU + r/fc2j), ki = const satisfying ki < k2 < ... < k^ 
(e.g., [38]). Taking n = 2, ki = 1.05, k2 = 1.2, and setting t = t — y, we obtain a 
two-soliton solution of the system (|5| shown in Fig. |2] 




Figure 2: Solution (15), (16) to the system ([s]) for v and w at t = 0. 

Further differential reductions of the system ^ can be obtained from the general 
approach of p9j. 
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3.3 Reduction to a pair of (l+l)-dimensional linearly degener- 
ate systems 

The method of hydro dynamic reductions, originally proposed in fTfi [TS] in the context of 
the dispersionless KP hierarchy, was subsequently generalized to a whole class of multi- 
dimensional quasilinear systems in [131 ttH IE] • We will use the example of system ^ to 
illustrate the method, see also [33|, HQ]. Let us seek A^-phase solutions of the system (|5| 
in the form 

v = v{R\...,R^), w = w{R\...,R^) (17) 

where the phases (Riemann invariants) R^{x,y,t) satisfy a pair of compatible (1 + 1)- 
dimensional systems 

(18) 



Rl - X'{R)Rl, R'y - fj'{R)Rl. 

.e : 



Recall that the commutativity conditions of the flows (18) take the form 



- A* /iJ - /i' 

7^ j, dj = d/djij, no summation. The substitution of the ansatz (17) into ^ implies 



(19) 



diW = —XdiV, /i* 



w 
A^ 



which, taking into account the commutativity condition (19), leads to dj log A' = dj logw, 
didjlogw = 0, didjV = 0, i j. These equations are straightforward to integrate: 

v = J2v'iR'), w = l[f{R'), X' = ^\R^)l[f{R'), 

k k k 

here f^{R^), r]^{R^) and ip^{R^) are arbitrary functions of one variable such that (?7*)' = 
— {f^y/if^'P^)- After a reparametrization i?* — )■ /'(-R*) this simplifies to 



w 



n< y = ^"mllR', (r/7 = -l/(i?V)- 



From now on we fix a particular choice ip^{R^) = 1/i?* under which both systems in (18) 
become linearly degenerate: 



Rl = (H R") Ri., ri = (J2r''U 

k^i kj^i kj^i 



(20) 



Once a solution to equations (20) is known, the variables v and w can be reconstructed 
by the formulae 

v = -J2r\ w = y[rk 

k k 

Thus, the (2 + 1) -dimensional system ^ is decoupled into a pair of (1 -|- l)-dimensional 
linearly degenerate systems (20) with an arbitrarily large number of dependent variables 
N. Since the gradient catastrophe for one-dimensional linearly degenerate hyperbolic 
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systems does not occur [SHI ESI [231 EZ] , generic solutions arising from this construction 
will not break down in finite time. The general solution of Eqs. (20) can be represented 
in the following parametric form, 



(21) 



R 



see [12], where ai'G ^ arbitrary functions of one variable. In the next subsections 
we utilize this representation to construct particular solutions of system ^ representing 
various multiple-wave interactions. 

3.3.1 Reduction to the Jacobi inversion problem 



With the choice = y P2N+iiO^ ^ = l;---!^; where P2n+i{0 is a polynomial of 



degree 2N + 1 in ^, Eqs. (21) reduce to the Jacobi inversion problem which can solved 
using multi-dimensional theta functions. 

The equation /x^ = ^ 11?=! ~ -^i) with Aj G M, z = 1, . . . , 2N, Aj 7^ A^ for i j and 
Aj 7^ defines a hyperelliptic Riemann surface C of genus N with simple branch points at 
the Ai, and at infinity. In addition we choose Aj = cj^ We order the branch points 
in a way that Ai < A2 < . . . < < . . . < \2n- The differentials 

Vk = , k = l,...,N, 

form a basis for the holomorphic one-forms on C. On this surface we introduce a basis of 
the homology, where the a-cycles are in the upper sheet each of them encircling exactly 
one of the cuts between [Ai,A2], [A3, A4],. . . ,[A2Ar-2, A2Ar_i]; all 6-cycles start at the cut 
[A27v,oo], and each of them intersects exactly one a-cycle such that Oj ■ bj = 6ij. A 
normalized basis of holomorphic differentials uji, i = 1, N — 1 is given by 

N 

uji = ^Cijr]j, i = l,...,N, 
3=1 

where the Cij are fixed by the normalization condition 

ujj = 2Tri6ij, i,j = l,...,N. 

^More general hyperelliptic Riemann surfaces can be used in this context, but since we only want to 
illustrate this approach, we took this choice leading to simple formulae. With this choice, v will be a 
solution to the KdV equation 



10 



With this normahzation the matrix 



Uj, i,j = 1,...,N, 

is a Riemann matrix, i.e., it is symmetric and has a negative definite real part (see for 
instance [2]). Therefore the theta function with half integer characteristics [a] = [a'; a"] G 
R"^^ defined by 

e[a](x|B) = exp f ^(n + a', B(n + a')) + (n + a', x + 27ria") j , 

where the angular brackets denote the Euclidean scalar product, is an entire function for 
all z e C^. 



With this setting the Jacobi inversion problem following from (21) can be written in 
the form 

Y ^i^^ ~ ^i{N-l)y + Cil(-l)^(y - t) =■ Zi. 

Notice that an arbitrary constant can be added to each Zi refiecting the freedom to 
choose a lower limit of the integrals on the left-hand side. Formulae for the symmetric 
functions Xl^Li ^^"^ Il^i ™^ terms of theta functions can be found for instance in 
[2]: the former appears in the hyperelliptic solutions to the KdV equation, the latter (as 
a logarithm) in the hyperelliptic solutions of the sine-Gordon equation, which is related 
to the sinh-Gordon equation by the transformation u — )■ iu. Thus not surprisingly both 
equations which appear as reductions of the studied systems also play a role in the context 
of hyperelliptic solutions. We have 

N N 

^7 = - J] i?' = -(c + In e[a](;2)), w = \[r' = d + d,thiQ[a]{z), (22) 

1=1 i=l 

where the arbitrary non-singular constant characteristic refiects the possibility to change 
z by an arbitrary constant, and where c and d are constants with respect to the physical 
coordinates, 

TV 



c = ^y^CiN f , d = d^t^nQ[a*] ( [ u] , 
2vn^ J a. /i \Jo J 

where [a*] is an arbitrary odd half integer characteristic. It can be seen that Vt + Wx = 



by construction. The validity of the second equation in (12) is checked by applying 
derivatives to the theta functions and computing these numerically. Since this relation 
is not implemented in the code, its validity provides a strong test on the quality of the 
numerics. The code typically satisfies this test to the order of machine precision. We use 
the code [TB] developed for hyperelliptic solutions to KdV, Kadomtsev-Petviashvili and 
other equations. In Fig. [3] one can see a typical genus 3 solution. As already mentioned, v 
solves the KdV equation in y and x, the time t playing just the role of a parameter here. 
In Fig. |4} the same situation as in Fig. |3]can be seen for t = 0.1 thus illustrating the time 
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Figure 3: Solutions v and w (22) for 
e, -3, -e, 0, 1, 1 + 2e, 4] for e = 1 at t = 
chosen to be = 1/2 and a'/ = for z 



0. 



= 3 on the surface with branch points [—3 — 
The characteristics of the theta function are 



1,2,3. 
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X 



Figure 4: Solutions v and w (22) for 
e, -3, -e, 0, 1, 1 + 2e, 4] for e = 1 at t = 
chosen to be = 1/2 and a'/ = for z 



3 on the surface with branch points [—3 — 
0.1. The characteristics of the theta function are 



1,2,3. 
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evolution of the pattern. If the e in Fig. [3] tends to 0, the surface degenerates which leads 
to the solitonic limit of the KdV solution v. This can be seen in Fig. [5] where a 3-soliton 
solution to KdV appears. The characteristic phase shift at the point of collision between 
the solitons can be clearly recognized. In this sense the solutions on a non-degenerate 





Figure 5: Solutions v and w (|22) for 

10-1^ at t 



3 on the surface with branch points [—3 — 
0. The characteristics of the theta function 



e, -3,-e,0, 1,1 + 2e,4] for e 
are chosen to be a[ = 1/2 and a'l = for i = 1, 2, 3. 

surface can be loosely interpreted as an infinite train of 3-solitons. Going to higher genus, 
one obtains A^-phase solutions which show with increasing N more and more structure, 
as can be seen for instance in Fig. [6] 

3.3.2 Parametric plots of 3-wave interactions 



In the case = 3 Eqs. (20) reduce to 

d1 d2 p3 d1 d2 d1 d3 pS p3 pi p2 p3 

JXj — it JX -ttx^ — it JX -tt^, iXj — It It it^ 

and 



(23) 



— {^B?-\-R^ — R^) R^., R^ — {^R^-\-R'^ — R^R'^) R^, R^ — {^R^-\-R^ — R^R^) R"^, (24) 
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Figure 6: Solutions v and w (22) for iV = 4 on the surface with branch points 
[—4, —3, —2, —1, 0, 1, 2, 3, 4] at t = 0. The characteristics of the theta function are chosen 
to be a[ = 1/2 and a'^ = for i = 1, 2, 3,4. 
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respectively. The formulae (21) simplify to 



„ R? „ ^ 



9(0 



Fi^ 



R2 ^3 



(25) 



where /, g, h are arbitrary functions. For our purposes, it will be more convenient to work 
with an equivalent parametric representation of the general solution of (23)-(24), namely, 

= F(a), = G{b), = H{c), 

where a, b, c are functions of x, y, t defined implicitly via 

^dc = -V, (26) 



J F{a)-l""^J G(fe)-l"'^^J H{c)-1' 



a + b + c 



here F, G, if are three arbitrary functions, and all integrals are understood as indefinite. 
The corresponding solution of ([s]) is given by the formulae 



-F{a)-G{b)-H{c), w = F{a)G{b)H{c). 



(27) 



Relations (26) can be derived from (25) by setting 



R^ = F{a), R^ = G{b), R^ = H{c) 



and 



(e - m 



i?2 



(e - m 

9i0 



respectively. Let us choose, for instance, three soliton-like profiles 



F{a) = a — 



cosh fxa 



G{b) = 13- 



cosh /i6 



Hie) 



1 - 



cosh /ic 



(28) 



where a, /3, 7 are distinct constants and e is a small parameter. In this case equations (26 ) 
take the form 



a 



13% , , 



a 



aa Bb 'yc , . 

+ 7;— ^ + — ^ + e(...) = -y, a + b + c=-t, 



7 



where dots denote certain explicit functions of a, 6, c which are bounded along with their 
first order derivatives. Thus, for sufficiently small e the change of variables a, 6, c — )■ x, y, t 



16 



is nonsingular and well-defined globally, so that the corresponding v and w will also be 
nonsingular. Expressing c in the form c = —t — a — b and substituting into the formulas 



(27) for V, w, and the first two equations (26), one obtains v,w,x,y as functions of a 



and b. Thus, one can construct parametric plots of v{x,y,t) and w{x,y,t). Figs. [7] 
a nd [8] show snapshots of v and w aX t = and t = 6 for the parameter values a = 1 — 

(l+tan((^+7r/3)(l+tan(0-7r/3)) ^ _ _ ^ (1+tan ^){l+ta,n{4)~TT/3)) ^ _ -|^_ / (1+tan 0)(l+tan(0+7i-/3)) 



1+tan <j) 



l+tan((/i+7r/3) 



l+tan((/i-7r/3) 



where = 7r/9, and /i = 1, e = 0.05. At t = initial profiles represent three solitons 
having a triple collision at the origin (parameters are chosen in such a way that solitons 
meet at equal angles vr/3). The triple collision then disintegrates into pairwise interactions, 
shown at t = 6. 




7 50 



Figure 7: Solution (27), (28) to the system (Is]) for t> at t = and t = 6. 



-0,4 V 



-0.51 



-0.6\^,-- 
-50 



^50 
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Figure 8: Solution (27), (28) to the system ([s]) for w at t = and t = 6. 



Unfortunately, for > 3 one cannot reduce Eqs. (21) to parametric formulas, and 
has to solve them numerically. 
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3.3.3 Plots of A^-wave interactions 



To obtain a more convenient form for explicit integration, we rewrite Eqs. (21) as 



/ F^-\p,)dp, + ■■■ + / F^-\pN)dp 
f F^{p^)dpi + ■•• + / F^^{pN)dpM = 0, k = 
Pi + ...+PN = {-lf{y-t). 



N 



(29) 



1,2, ...,iV 



here Fi{pi) are arbitrary functions, and all integrals are understood as indefinite. The 
corresponding solution of ^ is given by the formulae 



w 



W 



Formulae (21) and (29) are identified upon setting _R* = Fi{pi), f = pi. Equations 
(29) can be used to construct nonlinear interactions of soliton-like profiles by taking, 
e.g., Fi{pi) = + ^Qg^^l^.p. , where e is sufficiently small and the constants are chosen 
in such a way that equations (29) are uniquely solvable for pi{x,y,t). 

This approach can in principle be used to build hump-like initial data by letting 
tend to infinity, thus providing in a loose sense a non-linear generalization of Fourier 
analysis for the studied PDE. We illustrate this concept below by studying examples for 
finite N numerically. The idea is to construct for t = a central hump from N elementary 
solutions of the form of the KdV soliton. 

The angles between these elementary solutions are chosen to be equal. From the 
second relation in (20) it can be seen that a condition on these angles is equivalent to a 
condition on the a,. 



N 



N 



To ensure that w{x,y,0) < we choose Fi = ai + 
for i = N — 1, N. In addition we fix ^ 



tan ipi 



N 



cosh'' PiPi 



iV. (30) 



for i = 1, . . . , N — 2 and Fi = Ui 



and 11^=1 = ~1) i\ms assuring that 
we are studying solutions on the zero background for v and on the —1 background for w. 
The first N — 2 ai then follow from (30) via the quadratic relations, 

l/ai - ai 



idMipi, 



N 



The remaining two are a consequence of the condition of vanishing sum of the a, and 
product equal to —1, 



OiN~l,N 
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The ijji are chosen in a way to generate equal angles between the elementary solutions in 
the (x, I/)— plane for t = 0. For numerical reasons one has to avoid ip = 7c/2 since the 
tangent diverges there. To obtain finite ipi also for even N, we shift them by a common 
factor, thus 



- ) 

N odd, 

' N even. 



N-2 ' 2{N~2) 

With this choice we assure that the values for tanip are as small as possible. This is 
numerically important since the appear in an x Vandermonde determinant, i.e., 
there will be A^th powers of each aj. For the same reason we also have to choose the 
signs carefully to avoid too large or too small values of aAr_i jy. With the above choice 
of the ipi, i = 1, . . . , N — 2, we have assured that with some aj a solution to one of the 



equations (30), also —aj and ±aj are solutions to these equations (either for ipj or for 
vr — tpj). We fix the signs in a way to obtain a minimal Yli^'^ a negative ni=i^ ^i- 

This is done by alternating the signs in the a,, i = 1, . . . , Nmax {Nmax = N/2 — 1 for A^ 
even and N^ax = {N — l)/2 for A^ odd). The remaining signs are chosen in a way that 
the ttj corresponding to n — ipi is equal to Thus we assure that a^-i and are 

real and that l/laNmaxI ^ l"^*! ^ I^^Afmaa:!- We give the values of the for two examples 
below. Since aN^ax ~ 2(A^ — 2)/7i, we can in double precision reach values of A^ = 14 in 
which case the conditioning of the Vandermonde matrix of the is of the order of 10~^^ 
(the absolute value of the ratio of the smallest to the largest eigenvalue of the matrix) . To 
treat higher values of A^ more than double precision is needed which is, however, beyond 
the scope of this paper. 

The goal of this subsection is to construct for t = a single hump for v with value 
c = const at X = y = 0. To this end we choose e = c/{N — 2) since N — 2 elementary 
solutions will overlap at the origin of the (x, y) plane with the above choices. Thus we 
take e = c/{N — 2) and take c = 0.1. The constant c has to be chosen in a way that 
w remains negative which seems to be the case in the numerical experiments discussed 
below if c < 1. A smaller value of c leads to a more rapid convergence of the iteration to 
be described below. 



For the given form of the Fi, i = 1, A^, we have {n > 1) 

J FJ^{p,)dp, = dtp, + 1^ Q e^oll-n,{p,\ 
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N-2 ( N-1 

N-1 



where /„(p) = dp' / cosh^"p'. Equations (29) thus take the form (A; = 1, 2, . . . , — 3) 

N-2 I 2\ N-2-j J / \ 1 I N-2 I Ar-2 

«i Pi + ^- J h\Pi) I + «7v-im-i + pn = -y 

p^ + ...+p^ = {-lf{y-t), 




i=l 



This system can be solved iteratively. Let V be the Vandermonde determinant of the ttj. 
Then we solve iteratively for the vector p, where Pn+i follows from the nth iterate via 



Vpn+i = J{Pn) 



with 



etc, and 



N-2 N-1 . 
JliPn) = - E E ( J J ^^^t'^'^^i^Pi) + ^' 



1=1 j=l 



J(po) = (x,-2/,0,...,0,(-l)^(y-t)) 



The iteration is stopped once \pn+i — Pn\ < S (we use 6 = 10^^). For a value of c ^ 1, 
this iteration converges rapidly (typically less than 4 iterations). If c ~ 1 the iteration 
converges only slowly or not at all. 

The results of the above procedure for = 6 (four elementary solutions) for t = 
can be seen in Fig. [9) Since this approach is concerned with hump-like structures in v, 
the outcome for w is not obvious, but can be seen in the lower part. The central hump 



for V at the origin of the {x, y)-plane disintegrates with time as can be seen in Fig. 10 



Since we fixed the maximum of the absolute value of v at the center to be of the order 
of c, the relative maximum of the modulus of the elementary waves tends with large 
to 0. This is illustrated in Fig. 11 for = 13. The central hump for v at the origin of 
the (x, |/)-plane disintegrates with time as can be seen in Fig. 12 Compare this figure 
with the time evolution of hump-like initial data for v studied numerically in the next 
subsection. 



3.4 Cauchy problems and numerical evaluations 



We will use the Friedrichs-symmetrized form (10), which is well-posed for w < 0, for nu- 
merical simulations of system (|5]). The task is to investigate the conjecture that solutions 
to the system ([s]) for generic initial data do not develop a point of gradient catastrophe 
in finite time. 

To this end we consider initial data in the Schwartzian class of rapidly decreasing 
functions. Such functions can be treated for numerical purposes as effectively periodic if 
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Figure 9: Solutions v and w to the system (29) for = 6 at t = and jl = (1, 100, 1, 1). 
The values of fl have been chosen in a way that the lateral extension of the elementary 
waves is of a similar order. This implies a ~ (0.81, —2.77, —0.36, 1.23, —0.59, 1.69). 
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-50 -50 




-50 -50 



Figure 10: Solutions v and w to the system (29) for TV = 6 at t = 6 and jl = (1, 100, 1, 1). 



the period is chosen large enough that the studied functions decrease below the machine 
precision (which is of the order 10^^^ for double precision computations in MATLAB [29] ) 
at the ends of the intervals. The Gibbs phenomenon in such a case is of the order of the 
rounding error and thus negligible. This setting allows the use of Fourier spectral methods 
which are known for their excellent approximation properties for smooth functions. Thus 
we will approximate the spatial dependence of the functions v and w via discrete Fourier 
series. We find that the solutions to the system ([s]) stay in the Schwartzian class in 
contrast, e.g., to solutions to the dispersionless Kadomtsev-Petviashvili (dKP) equation, 
see pi] for a numerical study. 

The time dependence is treated via a standard fourth order Runge-Kutta scheme. The 
quality of the numerics is controlled via the conserved quantity 



(f ^ — (w — Woo))dxdy 



of system ([s]), where Woo w{x,y,t). Since the conservation of this quantity 

is not implemented in the code, it will not be exactly observed in actual numerical ex- 
periments due to unavoidable numerical errors. Thus the numerical conservation of / 
can serve as an indicator for the achieved accuracy. In our experiments the error in the 
relative conservation of this quantity was below 10^^. It is known (see for instance [22] ) 
that integral quantities overestimate the accuracy in a L^o sense of numerical solutions 
to a partial differential equations by one to two orders of magnitude. Thus the results 
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Figure 11: Solutions v and w to the system (29) for 
(100, 100, 10, 10^ 10, 10^°, 1, 10^ 10, 10^ 



13 at t = and pi 



100). The values of pi have been chosen in a way 
that the lateral extension of the elementary waves is of a similar order. This implies 



a 



:i.OO, -1.16, 0.73, -1.73, 0.39, -7.10, 



-0.14,2.58, 



-0.58,1.37,-0.86,-0.18,5.68) 
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presented below are accurate to at least 10~^ which is well beyond plotting accuracy. 
Comparing the behavior of the numerical solutions to solutions to the dKP equation in 
[21] which show blow up, we see in the experiments performed below no indication of an 
increasing gradient of the solutions. Thus it seems that at least for the class of initial 
data considered, the numerical solutions will have global existence in time. 

As a concrete example we study the initial data v{x,y,0) = exp(— (x^ + y^)) and 
w{x,y,0) = —1. It can be seen in Fig. 13 that the solution for v stays smooth and 
localized in both spatial directions. 




Figure 13: Solution v to the system ^ for the initial data v{x,y,0) = exp(— (x^ + y^)) 
and w{x,y,0) = —1 for several values of t. 



A similar behavior can be observed for the solution w in Fig. [14} This solution starts 
with w = —1 and develops peaks above and below this initial surface. It can be seen that 
both functions v and w approach an elliptic shape which could indicate the existence of 
an elliptic localized solution which acts as an attractor for a wide class of initial data. 

The elliptic shape seems to be also present if one starts with several humps. In Fig. [15] 
we consider initial data of the form v{x, y, 0) = exp(— ((x — 3)^+?/^))+exp(— ((x+3)^+y^)) 
and w{x,y,0) = —1. As can be seen both humps develop into an ellipse. The behavior 
for w is as follows: since it starts as identically constant, two humps of opposite sign are 
forming. Both develop then an elliptic shape as seen in Fig. [15] 
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t=0.0 



t=0.98 




y -10 -10 X y -10 -10 



Figure 14: Solution w to the system ^ for the initial data v{x,y,0) = exp(— (x^ + y^)) 
and w{x, y,0) = —1 for several values of t. 
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10 -10 y 

X 

Figure 15: Solutions v and w to the system ^ for the initial data v{x, y, 0) = exp(— ((x — 
3)2 + y^)) + exp(-((x + 3)2 + ?/2)) and w{x, y, 0) = -1 for t = 4. 
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4 Proof of Theorem 1 



In Sect. 4.1 we summarize the integrability conditions for two-component (2+l)-dimensional 
quasilinear systems. Their differential consequences are used in Sect. 4.2 to estabhsh The- 
orem 1. 



4.1 Integrability conditions for 2-component systems in 2 + 1 
dimensions 

Here we present the integrabihty conditions for 2-component quasihnear systems of the 
form ([2]). Setting u = (f , wY and diagonahzing the matrix A, which is always possible in 
the 2x2 case, one obtains the equations 




P q 



r s / \ w 



0. (31) 



y 



The integrability conditions constitute a complicated set of second order differential con- 
straints for the coefficients a, b, p, q, r, s which were derived in [H] based on the method 
of hydrodynamic reductions. 
Equations for a: 

_ qa^b^ + 2qal + (s - p)a^ay, - ra^ ^ a^r^ ^ 2a^p^ - a^p^, 



(a — b)q r q 

-), 

qciyby -|- (s p^Qj^b^ -\- TOjyj QjySyj (^wQu 



O-UI ~l" ,qw , '^W\ fnn\ 

av 7- + ay{ \ ), (32) 

a — b q r 



(a — b)r 

Equations for b: 



raj}yj + {p- s)a^by, + qb^ b 

wPv byTu 

(0 — a)q q r 

bvw = b j^l ^ + byj{— + — ), (33) 
b — a q r 

_ ra^b^ + 2rbl + jp - s)b^b^ - qbl b^^q^ 2b^s^ - b^s^ 

[b — a)r q r 
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Equations for p: 

Pvv — ^ , + + + 

(a — by r q 

^{2q^aw - 2a„qy, + a^jPy,) - b^py + 2ryay, - 2a^(s^, + Pv + r^) + 2=^(2p„a^ - a-^py,) 

b — a 

r,/ \ 0,vbw byjPy + C2Siu ~\~ Pw)(^v , I Qui , 

Pyyj = ^S-P)- — ^ — h (34) 

[a — b)'^ b — a \ q r J 



_ ^ qjawby - ttybw) + {s - p)ayjbyj {p - s)by,py - qb^p^ - 2rSu,aw - ra^Pw 

Pww — r ' \ ' 

[a — by [b — a)r 

Pv^w QwPw 

r q 
Equations for s: 



^ r{au,by - Gybu,) + (p - g)a^6^ {s - p)aySy, - ra^^Syj - 2qpyby - qb^s^ 



(a - by 



{a — b)q 



Z[p — S) -. — ^ - r Su, I I J , 

[a — b)^ a — b \q r 



(35) 



(a - by 



+ bySy) - a^Sy, + 2qyjby - 26^ (p^ + + ?„) + ^(2s 

a — b 
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Equations for q and r: 

(p - s)a^,b^ + q{a^b,aj - aM , qb^ + (s - p)b^ 

qrww + rqwin = ^[p - s) t-. hg r h 

[a — by r a — b 

I \ So-ui-^ui ~l~ 26^yP5jj -|- b^q^ (dw 26u;)g'u) 

s - p) + r + 

a — b a — b 

+ bv{2pu) + 2sw + Q'j;) — 2bw{rw + Pv + Sv) , 



a — b 



r 2 , Q 



q r 

qa^b^ + {s - p)a^byj + ra^b^ q^q^ p^Su, , 

= [s - p) + + + 

+ qpvbv 

r{a — b) 

ray^by, + {p - s)aybu] + qa^by r^ry, p^Sy, ^ 

Tvw = [P-S) \ \ h 

q{a — b)^ r q 

bwi'f'Qv + Qfv) + {p — s){aySy] + qPvbv 

q{b — a) ' 

{s - p)ayby + r{aybyj - ayjby) qy,ray, + {p - s)ay 

qrvv + rqvv = ^[s - p) — +r + 

[a — by q b — a 

/ \ 2byPy ~\~ ^dySy ~\~ dyVyj (by 2i(lyjry 

[p - s) ^ + q ^ + 

b — a b — a 

Jjyqy + ay,{2Sy + 2py + T y,) " 2 a^(g^ + Sy,+Py,) ^ 



(36) 



b — a 



fy + -Pvqw - rvQv + Pv{2sy + ryj); 
r q 

These formulas are completely symmetric under the identification v -H- w, a ^ b, p ^ 
s, q r. Our next goal is to rewrite these equations in a somewhat different form which 
will be used in the proof of Theorem 1. First of all, Eqs. (32)2 and (33)2 can be cast into 
the form 

d\nqr = n, ^] = f ^ + rft; + f ^ + ^f±^^ rfw;. (37) 

\by, a — b J \ dv b — a J 

In what follows we assume that all quantities ay, ay,, by, by, are nonzero, and parametrize 
the first order derivatives of q and r as 

Qw ^WW , byy by 

+ z n, — = h 



q ayj b — a r by a — b 

(3^ 

— = - — + T+m, — = + 7 +n, 

q by, by a — b r ay a^ b — a 



30 



where m and n are two auxiliary functions. Solving the remaining equations (32), (33) 
for the first order derivatives of p and s one gets 



Pv 



-qm 



a — b a — b 



(p-s), 



(^w , b^ byj 

— rn + q + [s - pj, 

tty b — a b — a 



Pv 



b — a 

by 
a — b 



[s-p) + 



(P - s ) + 



ay{a - bj 2 W 6„ 



1 



b^{b - a) 



q -\ — r In 



(39) 



+ 1 + 



m + 2 



duibv 1 . n 

1 H — \n + 2 



tty + by 

b — a 

O^W ~t~ byj 

a — b 



Substituting (39) into the equations for g^^, and one gets 



Inq^ 



Inr, 



1 —mn + 2 



(a - bf 



Qjybyj 



ay{a — b)q^ 



byw by,{ay + by) + ay,by ttyjby 6^^771 

1 —mn + 2- 



b — a 



(a - by 



bw{b - a)r' 



(40) 



respectively. Finally, differentiating the equation (39)i with respect to w and comparing 



the result with (|34|)2, as well as differentiating (39)2 with respect to v and comparing with 

by I Ctyu] O'W \ , (^W , o O'vby „ 

- —qm + 2nr- + 3g- — = 0, 

bw \ cLy a — b J b — a (a — oj^ 

(41) 



(35 )2, we obtain 
by 



-qm 



ttyj I b. 



Uy 



VW by \ , r\ , o ^wbw 



-rn rn + 2mq + 3r- 

ay J tty \by, b — J o " {tt — b)^ 



-'w 



0. 



Equations (37), (40) and (41) will be crucial for the proof of Theorem 1. 



4.2 Proof of Theorem 1 

For definiteness, let us consider a two-component (3 + l)-dimensional system of the form 
([!]). We assume that the generic matrix of the family aMj has a simple spectrum, and that 
the dispersion relation det(/i/ + aMj) = defines an irreducible quadric. Diagonalizing 
the first matrix, we obtain the equations 

ut + A{u)u^ + B{u)uy + B{u)u, = 

where u = {v, wY, and the matrices A, B, B have the form 



A 



a 



b 



B 



p q 



r s 



B 



p q 



r s 
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Due to the irreducibility of the dispersion relation, one cannot have g = g = Oorr = f = 0. 
Thus, by taking appropriate hnear combinations of B and B, one can assume q, r, q, r to 
be nonzero. Since all traveling wave reductions of our system must be integrable, the 
matrices B, B, and B + B must satisfy the (2 + l)-dimensional integrability conditions 
from Sect. 4.1 (with the same matrix A). There are three cases to consider: 
Case 1. A generic matrix of the family a A + /3B + 13 B has constant eigenvalues. In 
particular, all coefficients of the dispersion relation must be constant. This implies that 
a,b,p,s,p,s,qr,qf,qf + rq must be constant, and the integrability conditions of Sect. 
4.1, applied to pairs [A, B), {A, B) and {A, B + B), show that any such system must be 
linearizable by a transformation of the form v — )■ f{v), w — )■ g{w). 

Case 2. A generic matrix of the family aA + PB + PB is linearly degenerate. Then, 
due to the irreducibility of the dispersion relation, one cannot have a situation when a 
generic matrix of this family has one constant and one non-constant eigenvalue. Thus, we 
can assume that the matrix A is linearly degenerate, and has non-constant eigenvalues. 
In this case A can be brought to the form (M), and both matrices B and B must be of 
the form ([s]): see Example 1 of Sect. 2. Any such system is therefore a traveling wave 
reduction of the system ([5]). 

Case 3. A generic matrix of the family aA + PB + PB is genuinely nonlinear (again, due 
to the irreducibility of the dispersion relation, one cannot have a 'mixed' situation when 
a generic matrix of the family has one linearly degenerate, and one genuinely nonlinear 
eigenvalue). In particular, a„ 7^ and 7^ 0. In this case we will prove that matrices 
A, B and B must be linearly dependent, so that such system cannot be essentially (3 + 1)- 



dimensional. Applying the condition (37) to the matrices B, B and B + B one obtains 



dlnqr = dlnqr = dln{q + g)(r + r), 

so that q = Xq, r = fir where A and fi are constants. In what follows we assume that 7^ 
and by ^ (if both of these derivatives vanish, one can take a = v,b = w, however, with 
such A the integrability conditions from Sect. 4.1 are not compatible; similarly, due to the 
irreducibility of the dispersion relation one cannot have a situation when only one of these 



derivatives vanishes). If A = /x, the equations (39) imply p = Xp + ra + r], s = Xs + rb + r], 
so that B = XB + TA + riI , and our system becomes essentially (2 + l)-dimensional. Thus, 
we assume X fi. This means that, whenever one obtains a corollary of the integrability 
conditions which does not depend explicitly on p and s, and involves a, b, q, r only, it 
should still hold after the rescaling q — > Xq, r — )■ fir. In particular, rescaling in this way 



the relations (40) one obtains m = n = 0. This, however, is inconsistent with relations 



(41): the last term does not contain m,n, and cannot vanish. This finishes the proof. 



5 Concluding remarks 

Our results suggest that there must be a restriction on the dimensions in which integrable 
n-component systems can live. Namely, for any number of components n, there must be 
an integer d{n) such that there exist no n-component integrable systems of hydrodynamic 
type in the dimensions higher than d{n) + 1. We proved that d{2) = 5, moreover, there 
exists a unique two-component (5 -f- l)-dimensional integrable system. It would be of 
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interest to find a general formula for d{n), and classify the systems which attain an upper 
bound of this inequality. 
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